function metric(T, ut, ux, uy, utdt, uxdt, uydt, utdx, utdy, uxdx, uxdy, uydx, uydy, ...
  utdtdt, utdtdx, utdtdy, utdxdx, utdxdy, utdydy,...
  uxdtdt, uxdtdx, uxdtdy, uxdxdx, uxdxdy, uxdydy,...
  uydtdt, uydtdx, uydtdy, uydxdx, uydxdy, uydydy ...
  )
  % input -- T, ut, ux, uy at a single time step, of size GRID_SIZExGRID_SIZE


  global Dt DX DZ;
  global Z;
  GRID_SIZE = size(T, 1);
  Z_SIZE = size(Z,2);

  COMPUTE_REFERENCES = false;


  % g -- metric
  % gd -- metric derivatives
  % gam -- gamma
  % gamd -- gamma derivatives
  global g g_old g1;
  global detg gr;
  global gam gam_old;

  global gamd gd;

  global GUd Gdd;

  gd   = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4, 4);
  gamd = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4, 4, 4);
  Gdd = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4);

  % *************************************************************
  % *************************************************************
  %   precompute x-y quantities before the r loop
  % *************************************************************
  % *************************************************************
  if false
    ut =  ones(GRID_SIZE, GRID_SIZE);
    ux = zeros(GRID_SIZE, GRID_SIZE);
    uy = zeros(GRID_SIZE, GRID_SIZE);
    T = ones(GRID_SIZE, GRID_SIZE);
    utdt = zeros(GRID_SIZE, GRID_SIZE);
    utdx = zeros(GRID_SIZE, GRID_SIZE);
    utdy = zeros(GRID_SIZE, GRID_SIZE);
    uxdt = zeros(GRID_SIZE, GRID_SIZE);
    uxdx = zeros(GRID_SIZE, GRID_SIZE);
    uxdy = zeros(GRID_SIZE, GRID_SIZE);
    uydt = zeros(GRID_SIZE, GRID_SIZE);
    uydx = zeros(GRID_SIZE, GRID_SIZE);
    uydy = zeros(GRID_SIZE, GRID_SIZE);

    utdtdt = zeros(GRID_SIZE, GRID_SIZE);
    utdtdx = zeros(GRID_SIZE, GRID_SIZE);
    utdtdy = zeros(GRID_SIZE, GRID_SIZE);
    utdxdx = zeros(GRID_SIZE, GRID_SIZE);
    utdxdy = zeros(GRID_SIZE, GRID_SIZE);
    utdydy = zeros(GRID_SIZE, GRID_SIZE);

    uxdtdt = zeros(GRID_SIZE, GRID_SIZE);
    uxdtdx = zeros(GRID_SIZE, GRID_SIZE);
    uxdtdy = zeros(GRID_SIZE, GRID_SIZE);
    uxdxdx = zeros(GRID_SIZE, GRID_SIZE);
    uxdxdy = zeros(GRID_SIZE, GRID_SIZE);
    uxdydy = zeros(GRID_SIZE, GRID_SIZE);

    uydtdt = zeros(GRID_SIZE, GRID_SIZE);
    uydtdx = zeros(GRID_SIZE, GRID_SIZE);
    uydtdy = zeros(GRID_SIZE, GRID_SIZE);
    uydxdx = zeros(GRID_SIZE, GRID_SIZE);
    uydxdy = zeros(GRID_SIZE, GRID_SIZE);
    uydydy = zeros(GRID_SIZE, GRID_SIZE);
  end

  
  % *************************************************************
  % zeroth order
  % *************************************************************
  % u^mu u^nu
  utut = ut.*ut;
  utux = ut.*ux;
  utuy = ut.*uy;
  uxux = ux.*ux;
  uxuy = ux.*uy;
  uyuy = uy.*uy;

  % *************************************************************
  % first order
  % *************************************************************

  % partial_lam u^lam
  delu = utdt+uxdx+uydy;

  % partial_lam u^lam u_mu u_nu
  htt1 =  delu.*utut;
  htx1 = -delu.*utux; 
  hty1 = -delu.*utuy;
  hxx1 =  delu.*uxux;
  hxy1 =  delu.*uxuy;
  hyy1 =  delu.*uyuy;

  %          -         -        -        -
  htt2 = -2.*ut.*(ut.*utdt+ux.*utdx+uy.*utdy);

  hxx2 = -2.*ux.*(ut.*uxdt+ux.*uxdx+uy.*uxdy);
  hyy2 = -2.*uy.*(ut.*uydt+ux.*uydx+uy.*uydy);
  %      -           -        -        -     -
  htx2 =  (ux.*(ut.*utdt+ux.*utdx+uy.*utdy)+ut.*(ut.*uxdt+ux.*uxdx+uy.*uxdy));
  hty2 =  (uy.*(ut.*utdt+ux.*utdx+uy.*utdy)+ut.*(ut.*uydt+ux.*uydx+uy.*uydy));
  hxy2 = -(uy.*(ut.*uxdt+ux.*uxdx+uy.*uxdy)+ux.*(ut.*uydt+ux.*uydx+uy.*uydy));
  % without the -r up front
  htt2 = 2.*ut.*(ut.*utdt+utdx.*ux+utdy.*uy);
  htx2 = (-1).*ut.^2.*uxdt+(-1).*ux.*(utdx.*ux+utdy.*uy)+(-1).*ut.*(utdt.*ux+ux.*uxdx+uxdy.*uy);
  hty2 = (-1).*uy.*(utdx.*ux+utdy.*uy)+(-1).*ut.^2.*uydt+(-1).*ut.*(utdt.*uy+ux.*uydx+uy.*uydy);
  hxx2 = 2.*ux.*(ut.*uxdt+ux.*uxdx+uxdy.*uy);
  hxy2 = uxdy.*uy.^2+ut.*(uxdt.*uy+ux.*uydt)+ux.^2.*uydx+ux.*uy.*(uxdx+uydy);
  hyy2 = 2.*uy.*(ut.*uydt+ux.*uydx+uy.*uydy);


 
  % sigma_\mu_\nu

  stt = (1/2).*((-1).*utdt+2.*ut.^2.*utdt+(-1).*ut.^4.*utdt+ut.*utdx.*ux+(-1).*ut.^3.*utdx.*ux+(-1).*ut.*ux.*uxdt+ut.^3.*ux.*uxdt+uxdx+(-1).*ut.^2.*uxdx+ux.^2.*uxdx+ut.^2.*ux.^2.*uxdx+ut.*utdy.*uy+(-1).*ut.^3.*utdy.*uy+ux.*uxdy.*uy+ut.^2.*ux.*uxdy.*uy+(-1).*ut.*uy.*uydt+ut.^3.*uy.*uydt+ux.*uy.*uydx+ut.^2.*ux.*uy.*uydx+uydy+(-1).*ut.^2.*uydy+uy.^2.*uydy+ut.^2.*uy.^2.*uydy);
  stx = (1/2).*((-1).*utdx+ut.^2.*utdx+(-1).*ut.*utdt.*ux+ut.^3.*utdt.*ux+(-1).*utdx.*ux.^2+ut.^2.*utdx.*ux.^2+uxdt+(-1).*ut.^2.*uxdt+ux.^2.*uxdt+(-1).*ut.^2.*ux.^2.*uxdt+(-1).*ut.*ux.*uxdx+(-1).*ut.*ux.^3.*uxdx+(-1).*utdy.*ux.*uy+ut.^2.*utdy.*ux.*uy+(-1).*ut.*uxdy.*uy+(-1).*ut.*ux.^2.*uxdy.*uy+ux.*uy.*uydt+(-1).*ut.^2.*ux.*uy.*uydt+(-1).*ut.*uy.*uydx+(-1).*ut.*ux.^2.*uy.*uydx+ut.*ux.*uydy+(-1).*ut.*ux.*uy.^2.*uydy);
  sty = (1/2).*((-1).*utdy+ut.^2.*utdy+(-1).*ut.*ux.*uxdy+(-1).*ut.*utdt.*uy+ut.^3.*utdt.*uy+(-1).*utdx.*ux.*uy+ut.^2.*utdx.*ux.*uy+ux.*uxdt.*uy+(-1).*ut.^2.*ux.*uxdt.*uy+ut.*uxdx.*uy+(-1).*ut.*ux.^2.*uxdx.*uy+(-1).*utdy.*uy.^2+ut.^2.*utdy.*uy.^2+(-1).*ut.*ux.*uxdy.*uy.^2+uydt+(-1).*ut.^2.*uydt+uy.^2.*uydt+(-1).*ut.^2.*uy.^2.*uydt+(-1).*ut.*ux.*uydx+(-1).*ut.*ux.*uy.^2.*uydx+(-1).*ut.*uy.*uydy+(-1).*ut.*uy.^3.*uydy);

  sxx = (1/2).*((-1).*utdt+ut.^2.*utdt+(-1).*ut.*utdx.*ux+(-1).*utdt.*ux.^2+(-1).*ut.^2.*utdt.*ux.^2+(-1).*ut.*utdx.*ux.^3+ut.*ux.*uxdt+ut.*ux.^3.*uxdt+uxdx+2.*ux.^2.*uxdx+ux.^4.*uxdx+ut.*utdy.*uy+(-1).*ut.*utdy.*ux.^2.*uy+ux.*uxdy.*uy+ux.^3.*uxdy.*uy+(-1).*ut.*uy.*uydt+ut.*ux.^2.*uy.*uydt+ux.*uy.*uydx+ux.^3.*uy.*uydx+(-1).*uydy+(-1).*ux.^2.*uydy+(-1).*uy.^2.*uydy+ux.^2.*uy.^2.*uydy);
  sxy = (1/2).*((-1).*ut.*utdy.*ux+uxdy+ux.^2.*uxdy+(-1).*ut.*utdx.*uy+(-1).*utdt.*ux.*uy+(-1).*ut.^2.*utdt.*ux.*uy+(-1).*ut.*utdx.*ux.^2.*uy+ut.*uxdt.*uy+ut.*ux.^2.*uxdt.*uy+ux.*uxdx.*uy+ux.^3.*uxdx.*uy+(-1).*ut.*utdy.*ux.*uy.^2+uxdy.*uy.^2+ux.^2.*uxdy.*uy.^2+ut.*ux.*uydt+ut.*ux.*uy.^2.*uydt+uydx+ux.^2.*uydx+uy.^2.*uydx+ux.^2.*uy.^2.*uydx+ux.*uy.*uydy+ux.*uy.^3.*uydy);

  syy = (1/2).*((-1).*utdt+ut.^2.*utdt+ut.*utdx.*ux+(-1).*ut.*ux.*uxdt+(-1).*uxdx+(-1).*ux.^2.*uxdx+(-1).*ut.*utdy.*uy+ux.*uxdy.*uy+(-1).*utdt.*uy.^2+(-1).*ut.^2.*utdt.*uy.^2+(-1).*ut.*utdx.*ux.*uy.^2+ut.*ux.*uxdt.*uy.^2+(-1).*uxdx.*uy.^2+ux.^2.*uxdx.*uy.^2+(-1).*ut.*utdy.*uy.^3+ux.*uxdy.*uy.^3+ut.*uy.*uydt+ut.*uy.^3.*uydt+ux.*uy.*uydx+ux.*uy.^3.*uydx+uydy+2.*uy.^2.*uydy+uy.^4.*uydy);



  % *****************************************************
  %  second order
  % *****************************************************

  v3t = ut.*utdt.^2+(-1).*ut.^3.*utdt.^2+utdtdt+(-1).*ut.^2.*utdtdt+(-1).*ut.*utdtdx.*ux+utdt.*utdx.*ux+(-2).*ut.^2.*utdt.*utdx.*ux+(-1).*ut.*utdx.^2.*ux.^2+(-1).*ut.*uxdt.^2+ut.^3.*uxdt.^2+uxdtdx+(-1).*ut.^2.*uxdtdx+(-1).*ux.*uxdt.*uxdx+2.*ut.^2.*ux.*uxdt.*uxdx+ut.*ux.^2.*uxdx.^2+(-1).*ut.*ux.*uxdxdx+(-1).*ut.*utdtdy.*uy+utdt.*utdy.*uy+(-2).*ut.^2.*utdt.*utdy.*uy+(-2).*ut.*utdx.*utdy.*ux.*uy+(-1).*ut.*uxdxdy.*uy+(-1).*uxdt.*uxdy.*uy+2.*ut.^2.*uxdt.*uxdy.*uy+2.*ut.*ux.*uxdx.*uxdy.*uy+(-1).*ut.*utdy.^2.*uy.^2+ut.*uxdy.^2.*uy.^2+(-1).*ut.*uydt.^2+ut.^3.*uydt.^2+uydtdy+(-1).*ut.^2.*uydtdy+(-1).*ux.*uydt.*uydx+2.*ut.^2.*ux.*uydt.*uydx+ut.*ux.^2.*uydx.^2+(-1).*ut.*ux.*uydxdy+(-1).*uy.*uydt.*uydy+2.*ut.^2.*uy.*uydt.*uydy+2.*ut.*ux.*uy.*uydx.*uydy+ut.*uy.^2.*uydy.^2+(-1).*ut.*uy.*uydydy;
  v3x = utdtdx+ut.*utdt.*utdx+ut.^2.*utdt.^2.*ux+ut.*utdtdt.*ux+utdx.^2.*ux+utdtdx.*ux.^2+2.*ut.*utdt.*utdx.*ux.^2+utdx.^2.*ux.^3+(-1).*ut.^2.*ux.*uxdt.^2+ut.*ux.*uxdtdx+(-1).*ut.*uxdt.*uxdx+(-2).*ut.*ux.^2.*uxdt.*uxdx+(-1).*ux.*uxdx.^2+(-1).*ux.^3.*uxdx.^2+uxdxdx+ux.^2.*uxdxdx+utdx.*utdy.*uy+utdtdy.*ux.*uy+2.*ut.*utdt.*utdy.*ux.*uy+2.*utdx.*utdy.*ux.^2.*uy+ux.*uxdxdy.*uy+(-2).*ut.*ux.*uxdt.*uxdy.*uy+(-1).*uxdx.*uxdy.*uy+(-2).*ux.^2.*uxdx.*uxdy.*uy+utdy.^2.*ux.*uy.^2+(-1).*ux.*uxdy.^2.*uy.^2+(-1).*ut.^2.*ux.*uydt.^2+ut.*ux.*uydtdy+(-1).*ut.*uydt.*uydx+(-2).*ut.*ux.^2.*uydt.*uydx+(-1).*ux.*uydx.^2+(-1).*ux.^3.*uydx.^2+uydxdy+ux.^2.*uydxdy+(-2).*ut.*ux.*uy.*uydt.*uydy+(-1).*uy.*uydx.*uydy+(-2).*ux.^2.*uy.*uydx.*uydy+(-1).*ux.*uy.^2.*uydy.^2+ux.*uy.*uydydy;
  v3y = utdtdy+ut.*utdt.*utdy+utdx.*utdy.*ux+uxdxdy+(-1).*ut.*uxdt.*uxdy+(-1).*ux.*uxdx.*uxdy+ut.^2.*utdt.^2.*uy+ut.*utdtdt.*uy+utdy.^2.*uy+utdtdx.*ux.*uy+2.*ut.*utdt.*utdx.*ux.*uy+utdx.^2.*ux.^2.*uy+(-1).*ut.^2.*uxdt.^2.*uy+ut.*uxdtdx.*uy+(-2).*ut.*ux.*uxdt.*uxdx.*uy+(-1).*ux.^2.*uxdx.^2.*uy+ux.*uxdxdx.*uy+(-1).*uxdy.^2.*uy+utdtdy.*uy.^2+2.*ut.*utdt.*utdy.*uy.^2+2.*utdx.*utdy.*ux.*uy.^2+uxdxdy.*uy.^2+(-2).*ut.*uxdt.*uxdy.*uy.^2+(-2).*ux.*uxdx.*uxdy.*uy.^2+utdy.^2.*uy.^3+(-1).*uxdy.^2.*uy.^3+(-1).*ut.^2.*uy.*uydt.^2+ut.*uy.*uydtdy+(-2).*ut.*ux.*uy.*uydt.*uydx+(-1).*ux.^2.*uy.*uydx.^2+ux.*uy.*uydxdy+(-1).*ut.*uydt.*uydy+(-2).*ut.*uy.^2.*uydt.*uydy+(-1).*ux.*uydx.*uydy+(-2).*ux.*uy.^2.*uydx.*uydy+(-1).*uy.*uydy.^2+(-1).*uy.^3.*uydy.^2+uydydy+uy.^2.*uydydy;

  v4t = utdtdt+(-2).*ut.^2.*utdtdt+ut.^4.*utdtdt+(-2).*ut.*utdtdx.*ux+2.*ut.^3.*utdtdx.*ux+ut.^2.*utdxdx.*ux.^2+ut.*ux.*uxdtdt+(-1).*ut.^3.*ux.*uxdtdt+uxdtdx+(-1).*ut.^2.*uxdtdx+ux.^2.*uxdtdx+(-2).*ut.^2.*ux.^2.*uxdtdx+(-1).*ut.*ux.*uxdxdx+(-1).*ut.*ux.^3.*uxdxdx+(-2).*ut.*utdtdy.*uy+2.*ut.^3.*utdtdy.*uy+2.*ut.^2.*utdxdy.*ux.*uy+ux.*uxdtdy.*uy+(-2).*ut.^2.*ux.*uxdtdy.*uy+(-1).*ut.*uxdxdy.*uy+(-2).*ut.*ux.^2.*uxdxdy.*uy+ut.^2.*utdydy.*uy.^2+(-1).*ut.*ux.*uxdydy.*uy.^2+ut.*uy.*uydtdt+(-1).*ut.^3.*uy.*uydtdt+ux.*uy.*uydtdx+(-2).*ut.^2.*ux.*uy.*uydtdx+uydtdy+(-1).*ut.^2.*uydtdy+uy.^2.*uydtdy+(-2).*ut.^2.*uy.^2.*uydtdy+(-1).*ut.*ux.^2.*uy.*uydxdx+(-1).*ut.*ux.*uydxdy+(-2).*ut.*ux.*uy.^2.*uydxdy+(-1).*ut.*uy.*uydydy+(-1).*ut.*uy.^3.*uydydy;
  v4x = utdtdx+(-1).*ut.^2.*utdtdx+ut.*utdtdt.*ux+(-1).*ut.^3.*utdtdt.*ux+(-1).*ut.*utdxdx.*ux+utdtdx.*ux.^2+(-2).*ut.^2.*utdtdx.*ux.^2+(-1).*ut.*utdxdx.*ux.^3+ut.^2.*ux.^2.*uxdtdt+2.*ut.*ux.*uxdtdx+2.*ut.*ux.^3.*uxdtdx+uxdxdx+2.*ux.^2.*uxdxdx+ux.^4.*uxdxdx+(-1).*ut.*utdxdy.*uy+utdtdy.*ux.*uy+(-2).*ut.^2.*utdtdy.*ux.*uy+(-2).*ut.*utdxdy.*ux.^2.*uy+2.*ut.*ux.^2.*uxdtdy.*uy+2.*ux.*uxdxdy.*uy+2.*ux.^3.*uxdxdy.*uy+(-1).*ut.*utdydy.*ux.*uy.^2+ux.^2.*uxdydy.*uy.^2+ut.^2.*ux.*uy.*uydtdt+ut.*uy.*uydtdx+2.*ut.*ux.^2.*uy.*uydtdx+ut.*ux.*uydtdy+2.*ut.*ux.*uy.^2.*uydtdy+ux.*uy.*uydxdx+ux.^3.*uy.*uydxdx+uydxdy+ux.^2.*uydxdy+uy.^2.*uydxdy+2.*ux.^2.*uy.^2.*uydxdy+ux.*uy.*uydydy+ux.*uy.^3.*uydydy;
  v4y = utdtdy+(-1).*ut.^2.*utdtdy+(-1).*ut.*utdxdy.*ux+ut.*ux.*uxdtdy+uxdxdy+ux.^2.*uxdxdy+ut.*utdtdt.*uy+(-1).*ut.^3.*utdtdt.*uy+(-1).*ut.*utdydy.*uy+utdtdx.*ux.*uy+(-2).*ut.^2.*utdtdx.*ux.*uy+(-1).*ut.*utdxdx.*ux.^2.*uy+ut.^2.*ux.*uxdtdt.*uy+ut.*uxdtdx.*uy+2.*ut.*ux.^2.*uxdtdx.*uy+ux.*uxdxdx.*uy+ux.^3.*uxdxdx.*uy+ux.*uxdydy.*uy+utdtdy.*uy.^2+(-2).*ut.^2.*utdtdy.*uy.^2+(-2).*ut.*utdxdy.*ux.*uy.^2+2.*ut.*ux.*uxdtdy.*uy.^2+uxdxdy.*uy.^2+2.*ux.^2.*uxdxdy.*uy.^2+(-1).*ut.*utdydy.*uy.^3+ux.*uxdydy.*uy.^3+ut.^2.*uy.^2.*uydtdt+2.*ut.*ux.*uy.^2.*uydtdx+2.*ut.*uy.*uydtdy+2.*ut.*uy.^3.*uydtdy+ux.^2.*uy.^2.*uydxdx+2.*ux.*uy.*uydxdy+2.*ux.*uy.^3.*uydxdy+uydydy+2.*uy.^2.*uydydy+uy.^4.*uydydy;

  v5t = (-1).*(utdt+uxdx+uydy).*(ut.*utdt+(-1).*ut.^3.*utdt+utdx.*ux+(-1).*ut.^2.*utdx.*ux+ut.^2.*ux.*uxdt+ut.*ux.^2.*uxdx+utdy.*uy+(-1).*ut.^2.*utdy.*uy+ut.*ux.*uxdy.*uy+ut.^2.*uy.*uydt+ut.*ux.*uy.*uydx+ut.*uy.^2.*uydy);
  v5x = (utdt+uxdx+uydy).*((-1).*ut.^2.*utdt.*ux+(-1).*ut.*utdx.*ux.^2+ut.*uxdt+ut.*ux.^2.*uxdt+ux.*uxdx+ux.^3.*uxdx+(-1).*ut.*utdy.*ux.*uy+uxdy.*uy+ux.^2.*uxdy.*uy+ut.*ux.*uy.*uydt+ux.^2.*uy.*uydx+ux.*uy.^2.*uydy);
  v5y = (utdt+uxdx+uydy).*((-1).*ut.^2.*utdt.*uy+(-1).*ut.*utdx.*ux.*uy+ut.*ux.*uxdt.*uy+ux.^2.*uxdx.*uy+(-1).*ut.*utdy.*uy.^2+ux.*uxdy.*uy.^2+ut.*uydt+ut.*uy.^2.*uydt+ux.*uydx+ux.*uy.^2.*uydx+uy.*uydy+uy.^3.*uydy);

  v6t = (-1).*ut.*utdt.^2+ut.^3.*utdt.^2+(-1).*utdt.*utdx.*ux+2.*ut.^2.*utdt.*utdx.*ux+ut.*utdx.^2.*ux.^2+ut.*uxdt.^2+(-1).*ut.^3.*uxdt.^2+ux.*uxdt.*uxdx+(-2).*ut.^2.*ux.*uxdt.*uxdx+(-1).*ut.*ux.^2.*uxdx.^2+(-1).*utdt.*utdy.*uy+2.*ut.^2.*utdt.*utdy.*uy+2.*ut.*utdx.*utdy.*ux.*uy+uxdt.*uxdy.*uy+(-2).*ut.^2.*uxdt.*uxdy.*uy+(-2).*ut.*ux.*uxdx.*uxdy.*uy+ut.*utdy.^2.*uy.^2+(-1).*ut.*uxdy.^2.*uy.^2+ut.*uydt.^2+(-1).*ut.^3.*uydt.^2+ux.*uydt.*uydx+(-2).*ut.^2.*ux.*uydt.*uydx+(-1).*ut.*ux.^2.*uydx.^2+uy.*uydt.*uydy+(-2).*ut.^2.*uy.*uydt.*uydy+(-2).*ut.*ux.*uy.*uydx.*uydy+(-1).*ut.*uy.^2.*uydy.^2;
  v6x = (-1).*ut.*utdt.*utdx+(-1).*ut.^2.*utdt.^2.*ux+(-1).*utdx.^2.*ux+(-2).*ut.*utdt.*utdx.*ux.^2+(-1).*utdx.^2.*ux.^3+ut.^2.*ux.*uxdt.^2+ut.*uxdt.*uxdx+2.*ut.*ux.^2.*uxdt.*uxdx+ux.*uxdx.^2+ux.^3.*uxdx.^2+(-1).*utdx.*utdy.*uy+(-2).*ut.*utdt.*utdy.*ux.*uy+(-2).*utdx.*utdy.*ux.^2.*uy+2.*ut.*ux.*uxdt.*uxdy.*uy+uxdx.*uxdy.*uy+2.*ux.^2.*uxdx.*uxdy.*uy+(-1).*utdy.^2.*ux.*uy.^2+ux.*uxdy.^2.*uy.^2+ut.^2.*ux.*uydt.^2+ut.*uydt.*uydx+2.*ut.*ux.^2.*uydt.*uydx+ux.*uydx.^2+ux.^3.*uydx.^2+2.*ut.*ux.*uy.*uydt.*uydy+uy.*uydx.*uydy+2.*ux.^2.*uy.*uydx.*uydy+ux.*uy.^2.*uydy.^2;
  v6y = (-1).*ut.*utdt.*utdy+(-1).*utdx.*utdy.*ux+ut.*uxdt.*uxdy+ux.*uxdx.*uxdy+(-1).*ut.^2.*utdt.^2.*uy+(-1).*utdy.^2.*uy+(-2).*ut.*utdt.*utdx.*ux.*uy+(-1).*utdx.^2.*ux.^2.*uy+ut.^2.*uxdt.^2.*uy+2.*ut.*ux.*uxdt.*uxdx.*uy+ux.^2.*uxdx.^2.*uy+uxdy.^2.*uy+(-2).*ut.*utdt.*utdy.*uy.^2+(-2).*utdx.*utdy.*ux.*uy.^2+2.*ut.*uxdt.*uxdy.*uy.^2+2.*ux.*uxdx.*uxdy.*uy.^2+(-1).*utdy.^2.*uy.^3+uxdy.^2.*uy.^3+ut.^2.*uy.*uydt.^2+2.*ut.*ux.*uy.*uydt.*uydx+ux.^2.*uy.*uydx.^2+ut.*uydt.*uydy+2.*ut.*uy.^2.*uydt.*uydy+ux.*uydx.*uydy+2.*ux.*uy.^2.*uydx.*uydy+uy.*uydy.^2+uy.^3.*uydy.^2;



  
  for iz = 1:Z_SIZE
    z = Z(iz);
    % *************************************************************
    % zeroth order: 
    % *************************************************************
    % compute g/z^2 at zeroth order
    % compute drg/z^3
    f = (1-(z.*T./(1-z)).^3);
    fdz = (-3).*T.^3.*((-1)+z).^(-4).*z.^2;
    %fdr = 3.*(T./r).^3./r;

    % term numbers refer to eq (21) in Marks paper


    % term 1: -2 u_mu   dx^mu dr
    % term 1:  2 u_mu / z^2  dx^mu dr
    g(:,:,iz,1,4) = -ut; % gtr
    g(:,:,iz,2,4) =  ux; % gxr
    g(:,:,iz,3,4) =  uy; % gyr
    g(:,:,iz,4,4) =  0;

    % D[g/z^2,r] = -2 g/z^3 
    gd(:,:,iz,1,4,4) = -2.*g(:,:,iz,1,4);
    gd(:,:,iz,2,4,4) = -2.*g(:,:,iz,2,4);
    gd(:,:,iz,3,4,4) = -2.*g(:,:,iz,3,4);
    gd(:,:,iz,4,4,4) =  0;


    % term 2 only
    %  (1-z)^2/z^2  (1 - (z T/(1-z))^3)

    coefz   = -(1-z).^2.*f;
    coefzdz = -(1-z).^2.*fdz + 2.*(1-z).*f;

    g(:,:,iz,1,1) =  coefz.*utut;
    g(:,:,iz,1,2) = -coefz.*utux;
    g(:,:,iz,1,3) = -coefz.*utuy;
    g(:,:,iz,2,2) =  coefz.*uxux;
    g(:,:,iz,2,3) =  coefz.*uxuy;
    g(:,:,iz,3,3) =  coefz.*uyuy;

    gd(:,:,iz,1,1,4) = -2.*coefz.*utut + z.*coefzdz.*utut;
    gd(:,:,iz,1,2,4) =  2.*coefz.*utux - z.*coefzdz.*utux;
    gd(:,:,iz,1,3,4) =  2.*coefz.*utuy - z.*coefzdz.*utuy;
    gd(:,:,iz,2,2,4) = -2.*coefz.*uxux + z.*coefzdz.*uxux;
    gd(:,:,iz,2,3,4) = -2.*coefz.*uxuy + z.*coefzdz.*uxuy;
    gd(:,:,iz,3,3,4) = -2.*coefz.*uyuy + z.*coefzdz.*uyuy;

    % term 3 only
    % (1-z)^2/z^2
    coeffz   = (1-z).^2;
    coeffzdz = -2.*(1-z);
    g(:,:,iz,1,1) = g(:,:,iz,1,1) + coeffz.*( utut - 1);
    g(:,:,iz,1,2) = g(:,:,iz,1,2) + coeffz.*(-utux    );
    g(:,:,iz,1,3) = g(:,:,iz,1,3) + coeffz.*(-utuy    );
    g(:,:,iz,2,2) = g(:,:,iz,2,2) + coeffz.*( uxux + 1);
    g(:,:,iz,2,3) = g(:,:,iz,2,3) + coeffz.*( uxuy    );
    g(:,:,iz,3,3) = g(:,:,iz,3,3) + coeffz.*( uyuy + 1);

    gd(:,:,iz,1,1,4) = gd(:,:,iz,1,1,4) + (-2.*coeffz + z.*coeffzdz).*( utut - 1);
    gd(:,:,iz,1,2,4) = gd(:,:,iz,1,2,4) + (-2.*coeffz + z.*coeffzdz).*(-utux    );
    gd(:,:,iz,1,3,4) = gd(:,:,iz,1,3,4) + (-2.*coeffz + z.*coeffzdz).*(-utuy    );
    gd(:,:,iz,2,2,4) = gd(:,:,iz,2,2,4) + (-2.*coeffz + z.*coeffzdz).*( uxux + 1);
    gd(:,:,iz,2,3,4) = gd(:,:,iz,2,3,4) + (-2.*coeffz + z.*coeffzdz).*( uxuy    );
    gd(:,:,iz,3,3,4) = gd(:,:,iz,3,3,4) + (-2.*coeffz + z.*coeffzdz).*( uyuy + 1);





    % *************************************************************
    % first order: 
    % *************************************************************

    % remember g has 1/z^2 factored out!

    %  br = r./T;
    %  F,FF and Fdz,FFdz don't care about the 1/z^2 factor here yet
    F = (1/6).*(3.^(1/2).*pi+(-2).*3.^(1/2).*atan2(3.^(-1/2).*(2+((-2)+T).*z), z.*T)+3.*log(((-1)+z).^(-2).*(1+((-2)+T).*z+(1+(-1).*T+T.^2).*z.^2)));
    FF = (1-z).^2 ./ T.* F;

    Fdz = T.*((-1)+z).^(-1).*((-1)+z+(-1).*T.*z).*(1+((-2)+T).*z+(1+(-1).*T+T.^2).*z.^2).^(-1);
    FFdz = (-2.*(1-z) .* F + (1-z).^2 .* Fdz) ./ T;

    %  FF = 2.*r2./T.*F;

    %  Fdr = - (1 + br)./(br+br.^2+br.^3)./T;
    %  FFdr = 4.*r./T.*F + 2.*r2./T.*Fdr;

    %  (1-z)/z \partial_l u^l u_\mu u_nu
    %  = (1-z)z / z^2 
    % +356 
    if true
    g(:,:,iz,1,1) = g(:,:,iz,1,1) + (1-z).*z.*(htt1-htt2) + FF.*stt;
    g(:,:,iz,1,2) = g(:,:,iz,1,2) + (1-z).*z.*(htx1-htx2) + FF.*stx; 
    g(:,:,iz,1,3) = g(:,:,iz,1,3) + (1-z).*z.*(hty1-hty2) + FF.*sty;
    g(:,:,iz,2,2) = g(:,:,iz,2,2) + (1-z).*z.*(hxx1-hxx2) + FF.*sxx;
    g(:,:,iz,2,3) = g(:,:,iz,2,3) + (1-z).*z.*(hxy1-hxy2) + FF.*sxy;
    g(:,:,iz,3,3) = g(:,:,iz,3,3) + (1-z).*z.*(hyy1-hyy2) + FF.*syy;

    %  D[(1-z)/z] = -1/z^2 = -z/z^3
    %  % z derivatives
    gd(:,:,iz,1,1,4) = gd(:,:,iz,1,1,4) + -z.*(htt1-htt2) + (-2.*FF + z.*FFdz).*stt;
    gd(:,:,iz,1,2,4) = gd(:,:,iz,1,2,4) + -z.*(htx1-htx2) + (-2.*FF + z.*FFdz).*stx; 
    gd(:,:,iz,1,3,4) = gd(:,:,iz,1,3,4) + -z.*(hty1-hty2) + (-2.*FF + z.*FFdz).*sty;
    gd(:,:,iz,2,2,4) = gd(:,:,iz,2,2,4) + -z.*(hxx1-hxx2) + (-2.*FF + z.*FFdz).*sxx;
    gd(:,:,iz,2,3,4) = gd(:,:,iz,2,3,4) + -z.*(hxy1-hxy2) + (-2.*FF + z.*FFdz).*sxy;
    gd(:,:,iz,3,3,4) = gd(:,:,iz,3,3,4) + -z.*(hyy1-hyy2) + (-2.*FF + z.*FFdz).*syy;
    end



    % *************************************************************
    % second order: 
    % *************************************************************

    % - 2/br j(br) u_mu
    % j2v3 = -br./2;
    % j2v4 = (1/18).*(3.*br+(-6).*br.^2+2.*3.^(1/2).*pi.*br.^3+(-4).*3.^(1/2).*br.^3.*atan(3.^(-1/2).*(1+2.*br))+(sqrt(-1)*2).*3.^(1/2).*log(sqrt(-1)+(-1).*3.^(1/2)+(sqrt(-1)*2).*br)+(sqrt(-1)*(-2)).*3.^(1/2).*log(sqrt(-1)+3.^(1/2)+(sqrt(-1)*2).*br));
    % j2v5 = (1/18).*(3.*br+12.*br.^2+(-4).*3.^(1/2).*pi.*br.^3+8.*3.^(1/2).*br.^3.*atan(3.^(-1/2).*(1+2.*br))+(sqrt(-1)*(-4)).*3.^(1/2).*log(sqrt(-1)+(-1).*3.^(1/2)+(sqrt(-1)*2).*br)+(sqrt(-1)*4).*3.^(1/2).*log(sqrt(-1)+3.^(1/2)+(sqrt(-1)*2).*br));
    % j2v3 = real(j2v3);
    % j2v4 = real(j2v4);
    % j2v5 = real(j2v5);

    % j2v3dr = (-1/2);
    % j2v4dr = (1/2)+(-1).*r+3.^(-1/2).*pi.*r.^2+(-2).*3.^(-1/2).*r.^2.*atan(3.^(-1/2).*(1+2.*r));
    % j2v5dr = (-1/2)+2.*r+(-2).*3.^(-1/2).*pi.*r.^2+4.*3.^(-1/2).*r.^2.*atan(3.^(-1/2).*(1+2.*r));
    % 
    % % max(max(abs(imag(j2v5))))

    % j2t = j2v3 .* v3t  +  j2v4.*v4t  +  j2v5.*(v6t - v5t);
    % j2x = j2v3 .* v3x  +  j2v4.*v4x  +  j2v5.*(v6x - v5x);
    % j2y = j2v3 .* v3y  +  j2v4.*v4y  +  j2v5.*(v6y - v5y);

    % j2tdr = j2v3dr .* v3t  +  j2v4dr.*v4t  +  j2v5dr.*(v6t - v5t);
    % j2xdr = j2v3dr .* v3x  +  j2v4dr.*v4x  +  j2v5dr.*(v6x - v5x);
    % j2ydr = j2v3dr .* v3y  +  j2v4dr.*v4y  +  j2v5dr.*(v6y - v5y);

    % g(:,:,iz,1,1) = g(:,:,iz,1,1) - (-2) .* ut ./ br .* j2t;
    % g(:,:,iz,1,2) = g(:,:,iz,1,2) - (-2) .* ut ./ br .* j2x;
    % g(:,:,iz,1,3) = g(:,:,iz,1,3) - (-2) .* ut ./ br .* j2y;
    % g(:,:,iz,2,2) = g(:,:,iz,2,2) + (-2) .* ux ./ br .* j2x;
    % g(:,:,iz,2,3) = g(:,:,iz,2,3) + (-2) .* ux ./ br .* j2y;
    % g(:,:,iz,3,3) = g(:,:,iz,3,3) + (-2) .* uy ./ br .* j2y;

    % gd(:,:,iz,1,1,4) = gd(:,:,iz,1,1,4) - (-2) .* ut ./ br .* (j2tdr - j2t./r);
    % gd(:,:,iz,1,2,4) = gd(:,:,iz,1,2,4) - (-2) .* ut ./ br .* (j2xdr - j2x./r);
    % gd(:,:,iz,1,3,4) = gd(:,:,iz,1,3,4) - (-2) .* ut ./ br .* (j2ydr - j2y./r);
    % gd(:,:,iz,2,2,4) = gd(:,:,iz,2,2,4) + (-2) .* ux ./ br .* (j2xdr - j2x./r);
    % gd(:,:,iz,2,3,4) = gd(:,:,iz,2,3,4) + (-2) .* ux ./ br .* (j2ydr - j2y./r);
    % gd(:,:,iz,3,3,4) = gd(:,:,iz,3,3,4) + (-2) .* uy ./ br .* (j2ydr - j2y./r);

    % *************************************************************
    %  x & y derivatives -- all orders together 
    % *************************************************************
    % derivatives computed inside the r loop since T(x,y) is entangled in F(r/T)
    for i=1:4
      for j=i:4
        gd(:,:,iz,i,j,2) = dx2(reshape(g(:,:,iz,i,j), GRID_SIZE, GRID_SIZE));
        gd(:,:,iz,i,j,3) = dy2(reshape(g(:,:,iz,i,j), GRID_SIZE, GRID_SIZE));
      end
    end


  end % end of the z loop

  % fix the factors of 2 ------------ check with somebody TODO
  %   divide the upper triangular components of g (will symetrize below),
  %   and the upper traingular x,y,r derivaties (will also be symmetrized later)
  % if false
  %   for i=1:4
  %     for j=i+1:4
  %       g(:,:,:,i,j) = g(:,:,:,i,j)./2;
  %       gd(:,:,:,i,j,2) = gd(:,:,:,i,j,2)./2;
  %       gd(:,:,:,i,j,3) = gd(:,:,:,i,j,3)./2;
  %       gd(:,:,:,i,j,4) = gd(:,:,:,i,j,4)./2;
  %     end
  %   end
  % end

  % fill the symmetric components of g
  for i=1:4
    for j=i+1:4
      g(:,:,:,j,i) = g(:,:,:,i,j);
    end
  end
  
  % time derivative of g
  gd(:,:,:,:,:,1) = (g - g_old)./Dt;
  g_old = g;

  % pull of the overall 1/z^3 factor from the t,x,y derivatives!
  for iz = 1:Z_SIZE
    z = Z(iz);
    gd(:,:,iz,:,:,1) = gd(:,:,iz,:,:,1).*z;
    gd(:,:,iz,:,:,2) = gd(:,:,iz,:,:,2).*z;
    gd(:,:,iz,:,:,3) = gd(:,:,iz,:,:,3).*z;
  end

  % all derivatives of g computed at this time. fill the symm part
  for i=1:4
    for j=i+1:4
      gd(:,:,:,j,i,:) = gd(:,:,:,i,j,:);
    end
  end

  if (COMPUTE_REFERENCES)
    gref = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4);
    giref = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4);
    gdref = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4,4,4);
    for iz = 1:Z_SIZE
      z = Z(iz);
      gref(:,:,iz,1,1) = (-1).*((-1)+z).^(-1).*(((-1)+z).^3+T.^3.*z.^3);
      gref(:,:,iz,1,4) = -1;
      gref(:,:,iz,4,1) = -1;
      gref(:,:,iz,2,2) = ((-1)+z).^2;
      gref(:,:,iz,3,3) = ((-1)+z).^2;
      
      gdref(:,:,iz,1,1,4) = ((-1)+z).^(-2).*(2+(-6).*z+6.*z.^2+((-2)+T.^3).*z.^3);
      gdref(:,:,iz,1,4,4) = 2;
      gdref(:,:,iz,4,1,4) = 2;
      gdref(:,:,iz,2,2,4) = 2.*((-1)+z);
      gdref(:,:,iz,3,3,4) = 2.*((-1)+z);

      % giref(:,:,iz, 1,4) = 1;
      % giref(:,:,iz, 4,1) = 1;
      % giref(:,:,iz, 2,2) = r.^(-2);
      % giref(:,:,iz, 3,3) = r.^(-2);
      % giref(:,:,iz, 4,4) = r.^(-1).*(r+(-1).*T).*(r.^2+r.*T+T.^2);
    end

    gerr  = sum(sum(sum(sum(sum(sum(  (g-gref ).^2 ))))));
    gderr = sum(sum(sum(sum(sum(sum(  (gd-gdref ).^2 ))))));
  end



  

  % inverse of g, computed as z^2 gi
  gi = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4);

  detg = g(:,:,:,1,3).^2.*g(:,:,:,2,4).^2+(-2).*g(:,:,:,1,4).*g(:,:,:,1,3).*g(:,:,:,2,4).*g(:,:,:,2,3)+g(:,:,:,1,4).^2.*g(:,:,:,2,3).^2+(-2).*g(:,:,:,1,2).*g(:,:,:,1,3).*g(:,:,:,2,4).*g(:,:,:,3,4)+2.*g(:,:,:,1,4).*g(:,:,:,1,3).*g(:,:,:,2,2).*g(:,:,:,3,4)+(-2).*g(:,:,:,1,4).*g(:,:,:,1,2).*g(:,:,:,2,3).*g(:,:,:,3,4)+2.*g(:,:,:,1,1).*g(:,:,:,2,4).*g(:,:,:,2,3).*g(:,:,:,3,4)+g(:,:,:,1,2).^2.*g(:,:,:,3,4).^2+(-1).*g(:,:,:,1,1).*g(:,:,:,2,2).*g(:,:,:,3,4).^2+2.*g(:,:,:,1,4).*g(:,:,:,1,2).*g(:,:,:,2,4).*g(:,:,:,3,3)+(-1).*g(:,:,:,1,1).*g(:,:,:,2,4).^2.*g(:,:,:,3,3)+(-1).*g(:,:,:,1,4).^2.*g(:,:,:,2,2).*g(:,:,:,3,3);

  gi(:,:,:,1,1) = 2.*g(:,:,:,2,4).*g(:,:,:,2,3).*g(:,:,:,3,4)+(-1).*g(:,:,:,2,2).*g(:,:,:,3,4).^2+(-1).*g(:,:,:,2,4).^2.*g(:,:,:,3,3);
  gi(:,:,:,1,2) = (-1).*g(:,:,:,1,3).*g(:,:,:,2,4).*g(:,:,:,3,4)+(-1).*g(:,:,:,1,4).*g(:,:,:,2,3).*g(:,:,:,3,4)+g(:,:,:,1,2).*g(:,:,:,3,4).^2+g(:,:,:,1,4).*g(:,:,:,2,4).*g(:,:,:,3,3);
  gi(:,:,:,1,3) = g(:,:,:,1,3).*g(:,:,:,2,4).^2+(-1).*g(:,:,:,1,4).*g(:,:,:,2,4).*g(:,:,:,2,3)+(-1).*g(:,:,:,1,2).*g(:,:,:,2,4).*g(:,:,:,3,4)+g(:,:,:,1,4).*g(:,:,:,2,2).*g(:,:,:,3,4);
  gi(:,:,:,1,4) = (-1).*g(:,:,:,1,3).*g(:,:,:,2,4).*g(:,:,:,2,3)+g(:,:,:,1,4).*g(:,:,:,2,3).^2+g(:,:,:,1,3).*g(:,:,:,2,2).*g(:,:,:,3,4)+(-1).*g(:,:,:,1,2).*g(:,:,:,2,3).*g(:,:,:,3,4)+g(:,:,:,1,2).*g(:,:,:,2,4).*g(:,:,:,3,3)+(-1).*g(:,:,:,1,4).*g(:,:,:,2,2).*g(:,:,:,3,3);
  gi(:,:,:,2,2) = 2.*g(:,:,:,1,4).*g(:,:,:,1,3).*g(:,:,:,3,4)+(-1).*g(:,:,:,1,1).*g(:,:,:,3,4).^2+(-1).*g(:,:,:,1,4).^2.*g(:,:,:,3,3);
  gi(:,:,:,2,3) = (-1).*g(:,:,:,1,4).*g(:,:,:,1,3).*g(:,:,:,2,4)+g(:,:,:,1,4).^2.*g(:,:,:,2,3)+(-1).*g(:,:,:,1,4).*g(:,:,:,1,2).*g(:,:,:,3,4)+g(:,:,:,1,1).*g(:,:,:,2,4).*g(:,:,:,3,4);
  gi(:,:,:,2,4) = g(:,:,:,1,3).^2.*g(:,:,:,2,4)+(-1).*g(:,:,:,1,4).*g(:,:,:,1,3).*g(:,:,:,2,3)+(-1).*g(:,:,:,1,2).*g(:,:,:,1,3).*g(:,:,:,3,4)+g(:,:,:,1,1).*g(:,:,:,2,3).*g(:,:,:,3,4)+g(:,:,:,1,4).*g(:,:,:,1,2).*g(:,:,:,3,3)+(-1).*g(:,:,:,1,1).*g(:,:,:,2,4).*g(:,:,:,3,3);
  gi(:,:,:,3,3) = 2.*g(:,:,:,1,4).*g(:,:,:,1,2).*g(:,:,:,2,4)+(-1).*g(:,:,:,1,1).*g(:,:,:,2,4).^2+(-1).*g(:,:,:,1,4).^2.*g(:,:,:,2,2);
  gi(:,:,:,3,4) = (-1).*g(:,:,:,1,2).*g(:,:,:,1,3).*g(:,:,:,2,4)+g(:,:,:,1,4).*g(:,:,:,1,3).*g(:,:,:,2,2)+(-1).*g(:,:,:,1,4).*g(:,:,:,1,2).*g(:,:,:,2,3)+g(:,:,:,1,1).*g(:,:,:,2,4).*g(:,:,:,2,3)+g(:,:,:,1,2).^2.*g(:,:,:,3,4)+(-1).*g(:,:,:,1,1).*g(:,:,:,2,2).*g(:,:,:,3,4);
  gi(:,:,:,4,4) = (-1).*g(:,:,:,1,3).^2.*g(:,:,:,2,2)+2.*g(:,:,:,1,2).*g(:,:,:,1,3).*g(:,:,:,2,3)+(-1).*g(:,:,:,1,1).*g(:,:,:,2,3).^2+(-1).*g(:,:,:,1,2).^2.*g(:,:,:,3,3)+g(:,:,:,1,1).*g(:,:,:,2,2).*g(:,:,:,3,3);

  % divide by the denominator
  for i=1:4
    for j=i:4
      gi(:,:,:,i,j) = gi(:,:,:,i,j)./detg;
    end
  end

  % fill the symm part
  for i=1:4
    for j=i+1:4
      gi(:,:,:,j,i) = gi(:,:,:,i,j);
    end
  end

  %gierr = sum(sum(sum(sum(sum(sum(  (gi-giref ).^2 ))))))
  %gi = giref;


  % compute g*gi to check if the inverse works
  %   checked -- ok
  % dij = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4);
  % err = 0;
  % for i=1:4
  %   for j=1:4
  %     tmp = 0;
  %     for k=1:4
  %       tmp = tmp + gi(:,:,:,i,k).*g(:,:,:,k,j);
  %     end
  %     dij(:,:,:,i,j) = tmp; 
  %     if i==j
  %       err = err + (tmp-1).^2;
  %     else
  %       err = err + tmp.^2;
  %     end
  %   end
  % end
  % max(max(max(sqrt(err))))
  %

  

  
  % compute Gamma's
  % compute the symmetric part twice, what the hell, it's just addition
  %  NO LONGER TRUE: be really careful with the contractions: gi is z^2, g is 1/z^2, but gd has 1/z^3 only in the dz component :) 
  %  gd - 1/z^3
  %  gi - z^2
  % hence, gam -- 1/z
  for i=1:4
    for j=1:4
      for k=1:4
        gam(:,:,:,i,j,k) = 1/2*( ...
        gi(:,:,:,i,1).*(gd(:,:,:,1,j,k) + gd(:,:,:,1,k,j) - gd(:,:,:,j,k,1)) + ...
        gi(:,:,:,i,2).*(gd(:,:,:,2,j,k) + gd(:,:,:,2,k,j) - gd(:,:,:,j,k,2)) + ...
        gi(:,:,:,i,3).*(gd(:,:,:,3,j,k) + gd(:,:,:,3,k,j) - gd(:,:,:,j,k,3)) + ...
        gi(:,:,:,i,4).*(gd(:,:,:,4,j,k) + gd(:,:,:,4,k,j) - gd(:,:,:,j,k,4)) );
      end
    end
  end


  if (COMPUTE_REFERENCES)
    gamref  = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4,4,4);
    gamdref = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4,4,4,4);
    for iz = 1:Z_SIZE
      z = Z(iz);

      gamref(:,:,iz,1,1,1) = (1/2).*((-1)+z).^(-3).*(2.*(1+(-1).*z).*((-1)+z).^3+T.^3.*((-1)+z).*z.^3);
      gamref(:,:,iz,1,2,2) = (-1)+z;
      gamref(:,:,iz,1,3,3) = (-1)+z;
      gamref(:,:,iz,2,2,4) = ((-1)+z).^(-1);
      gamref(:,:,iz,2,4,2) = ((-1)+z).^(-1);
      gamref(:,:,iz,3,3,4) = ((-1)+z).^(-1);
      gamref(:,:,iz,3,4,3) = ((-1)+z).^(-1);
      gamref(:,:,iz,4,1,1) = (-1/2).*((-1)+z).^(-3).*(2.*(1+(-1).*z).*((-1)+z).^5+(-1).*T.^3.*((-1)+z).^3.*z.^3+T.^6.*z.^6);
      gamref(:,:,iz,4,1,4) = (-1/2).*((-1)+z).^(-2).*((-2).*((-1)+z).^3+T.^3.*z.^3);
      gamref(:,:,iz,4,4,1) = (-1/2).*((-1)+z).^(-2).*((-2).*((-1)+z).^3+T.^3.*z.^3);
      gamref(:,:,iz,4,2,2) = (1/2).*((-1)+z).^(-3).*(2.*(1+(-1).*z).*((-1)+z).^5+(-2).*T.^3.*((-1)+z).^3.*z.^3);
      gamref(:,:,iz,4,3,3) = (1/2).*((-1)+z).^(-3).*((-2)+12.*z+(-30).*z.^2+40.*z.^3+(-2).*T.^3.*((-1)+z).^3.*z.^3+(-30).*z.^4+12.*z.^5+(-2).*z.^6);
      gamref(:,:,iz,4,4,4) = -2;

      gamdref(:,:,iz,1,1,1,4) = (-1).*((-1)+z).^(-4).*(((-1)+z).^4+T.^3.*((-1)+z).*z.^3);
      gamdref(:,:,iz,1,2,2,4) = 1;
      gamdref(:,:,iz,1,3,3,4) = 1;
      gamdref(:,:,iz,2,2,4,4) = ((-1)+z).^(-3).*((-1)+3.*z+(-2).*z.^2);
      gamdref(:,:,iz,2,4,2,4) = ((-1)+z).^(-3).*((-1)+3.*z+(-2).*z.^2);
      gamdref(:,:,iz,3,3,4,4) = ((-1)+z).^(-3).*((-1)+3.*z+(-2).*z.^2);
      gamdref(:,:,iz,3,4,3,4) = ((-1)+z).^(-3).*((-1)+3.*z+(-2).*z.^2);
      gamdref(:,:,iz,4,1,1,4) = (1/2).*((-1)+z).^(-4).*(2.*T.^3.*((-1)+z).^4.*z.^3+T.^6.*(5+(-2).*z).*z.^6+(-2).*((-1)+z).^5.*(1+z+(-2).*z.^2));
      gamdref(:,:,iz,4,1,4,4) = ((-1)+z).^(-3).*((-1)+3.*z+(-3).*z.^2+z.^3+T.^3.*z.^3);
      gamdref(:,:,iz,4,4,1,4) = ((-1)+z).^(-3).*((-1)+3.*z+(-3).*z.^2+z.^3+T.^3.*z.^3);
      gamdref(:,:,iz,4,2,2,4) = (1/2).*((-1)+z).^(-4).*((-4).*T.^3.*((-1)+z).^4.*z.^3+2.*((-1)+z).^5.*(1+z+(-2).*z.^2));
      gamdref(:,:,iz,4,3,3,4) = (1/2).*((-1)+z).^(-4).*((-2)+8.*z+(-6).*z.^2+(-20).*z.^3+(-4).*T.^3.*((-1)+z).^4.*z.^3+50.*z.^4+(-48).*z.^5+22.*z.^6+(-4).*z.^7);
      gamdref(:,:,iz,4,4,4,4) = 2;
    end

    %gamerr = sum(sum(sum(sum(sum(sum( (gam-gamref).^2 ))))))
    %gam = gamref;
  end

  % compute Gamma's derivatives..... errrr^2

  % time derivative
  gamd(:,:,:,:,:,:,1) = (gam - gam_old) ./Dt;
  gam_old = gam;

  % space derivative
  %  for loops only through 1/2 of the symmetric part
  for iz=1:Z_SIZE
    for i=1:4
      for j=1:4
        for k=1:4 % to compute one part of the symmetric tensor change 1 to j here, and uncomment the portion below
          gamd(:,:,iz,i,j,k,2) = dx2(reshape(gam(:,:,iz,i,j,k), GRID_SIZE, GRID_SIZE));
          gamd(:,:,iz,i,j,k,3) = dy2(reshape(gam(:,:,iz,i,j,k), GRID_SIZE, GRID_SIZE));
        end
      end
    end
  end
  % fill the symmetric part
  % for i=1:4
  %   for j=1:4
  %     for k=j+1:4
  %       gamd(:,:,:,i,k,j,2) = gamd(:,:,:,i,j,k,2); ;
  %       gamd(:,:,:,i,k,j,3) = gamd(:,:,:,i,j,k,3); ;
  %     end
  %   end
  % end

  % pull out an overall 1/z^2 from the t,x,y derivatives:
  for iz=1:Z_SIZE
    z = Z(iz);
    gamd(:,:,iz,:,:,:,1) = gamd(:,:,iz,:,:,:,1).*z;
    gamd(:,:,iz,:,:,:,2) = gamd(:,:,iz,:,:,:,2).*z;
    gamd(:,:,iz,:,:,:,3) = gamd(:,:,iz,:,:,:,3).*z;
  end

  % compute the r derivative, keeping in mind that gamma has 1/z factored out. Thus gamdz is 1/z^2
  for iz=2:Z_SIZE
    z=Z(iz);
    gamd(:,:,iz,:,:,:,4) = -gam(:,:,iz,:,:,:) + z.* ( gam(:,:,iz,:,:,:) - gam(:,:,iz-1,:,:,:) ) ./ ( Z(iz) - Z(iz-1) ) ;
  end
  %gamd(:,:,1,:,:,:,4) = gamdref(:,:,1,:,:,:,4);

  if (COMPUTE_REFERENCES)
    gamderr = sum(sum(sum(sum(sum(sum(sum(  (gamd-gamdref).^2 )))))))*DZ
  end

  % hurray. 


  % compute the (drumroll) RRRic-ci Scalarrrr!! (applause)

  gr = 0;
  for a=1:4
    for b=1:4
      for c=1:4
        tmp = 0;
        for d=1:4
          tmp = tmp + gam(:,:,:,d,a,b).*gam(:,:,:,c,c,d) - gam(:,:,:,d,a,c).*gam(:,:,:,c,b,d);
        end
        gr = gr + gi(:,:,:,a,b).*( gamd(:,:,:,c,a,b,c) - gamd(:,:,:,c,a,c,b) + tmp );
      end
    end
  end

  gr(:,:,1) = -12;
  %meangr = sum(sum(sum(gr)))/GRID_SIZE/GRID_SIZE/Z_SIZE


  % next up: RUddd
  %  ruddd = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4, 4, 4);
  %  for a=1:4
  %    for b=1:4
  %      for c=1:4
  %        for d=1:4
  %          tmp = 0;
  %          for k=1:4
  %            tmp = tmp + gam(:,:,:,a,c,k) .* gam(:,:,:,k,d,b) - gam(:,:,:,a,d,k) .* gam(:,:,:,k,c,b);
  %          end
  %          ruddd(:,:,:,a,b,c,d) = gamd(:,:,:,a,d,b,c) - gamd(:,:,:,a,c,b,d) + tmp;
  %        end
  %      end
  %    end
  %  end

  %  % Rdd
  %  rdd = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4);
  %  for k=1:4
  %    rdd = rdd + reshape(ruddd(:,:,:,k,:,k,:), GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4);
  %  end



  %  gr2 = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE);
  %  for i=1:4
  %    for j=1:4
  %      gr2 = gr2 + gi(:,:,:,i,j).*rdd(:,:,:,j,i);
  %    end
  %  end

  %  gr = gr2;

  %  Gdd = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4);
  %  for i=1:4
  %    for j=1:4
  %      Gdd(:,:,:,i,j) = rdd(:,:,:,i,j) - (1./2 .* gr + 3).* g(:,:,:,i,j);
  %    end
  %  end

  %gr = gd(:,:,:,3,3,4);



  function r = dx2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for y = 1:GRID_SIZE
      r(1,y) = - arg(GRID_SIZE,y)/2 + arg(2, y)/2;
    end
    for y = 1:GRID_SIZE
      r(GRID_SIZE,y) =  -arg(GRID_SIZE-1,y)/2 + arg(1, y)/2;
    end
    for x = 2:GRID_SIZE-1
      for y = 1:GRID_SIZE
        r(x, y) = (arg(x+1,y) - arg(x-1,y))/2;   
      end
    end 
    r = r / DX;
  end
  function r = dy2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for x = 1:GRID_SIZE
      r(x,1) = - arg(x,GRID_SIZE)/2 + arg(x,2)/2;
    end
    for x = 1:GRID_SIZE
      r(x,GRID_SIZE) =  -arg(x,GRID_SIZE-1)/2 + arg(x,1)/2;
    end
    for x = 1:GRID_SIZE
      for y = 2:GRID_SIZE-1
        r(x, y) = (arg(x,y+1) - arg(x,y-1))/2;   
      end
    end 
    r = r / DX;
  end


end

  function r = dx(arg)
    global d1;
    r = d1 * arg;
  end
  function r = dy(arg)
    global d1;
    r = arg * (d1.');
  end
  function r = dxdx(arg)
    global d2;
    r = d2 * arg;
  end
  function r = dydy(arg)
    global d2;
    r = arg * (d2.');
  end
